/*
Date: October 2025
Project: Income and Child Maltreatment: Evidence from a Discontinuity in Tax Benefits
Author: Katherine Rittenhouse
Purpose: This file conducts analyses for various figures and tables throughout the paper and supplementary appendices
Files in: analysis.dta
*/

clear all
set more off

set maxvar 32767


*panelcombine program
include "https://raw.githubusercontent.com/steveofconnell/PanelCombine/master/PanelCombine.do"

log using "analysis"
use "analysis.dta" if abs(et)<=60,clear
local bw 60
local lower -`bw'
local upper `bw'
local tick 10


local donutl 8
local donutr 8
local fmt tex

replace et = et*-1

*instrument = born in december
gen instrument = (et>0)
gen et_r = et*instrument
gen transformed = et/(`bw')
gen inbandwidth = abs(et)<`bw'
gen triwt = (1-abs(transformed))*inbandwidth
gen indonut = 0

replace indonut = 1 if et<-(`donutl'-1) | et>(`donutr')

keep if recentered<=2017

gen patmiss = 0
replace patmiss = 1 if feduc==9
replace patmiss = 1 if dadrace==9
replace patmiss = 1 if dadagebin==9

**BALANCE**

gen Mom_Black = momrace==2
gen Mom_White = momrace==1
gen Mom_RaceOther = (momrace==3 | momrace==4)
gen Mom_Hispanic = momrace==5
gen Mom_RaceMiss = momrace==9
gen Mom_LowEd = meduc==1
gen Mom_HighEd = meduc==2
gen Mom_EdMiss = meduc==9
gen Mom_Age = momage
gen Mom_AgeMiss = momagebin==9
gen Mom_AnyMissing = max(Mom_AgeMiss,Mom_RaceMiss,Mom_EdMiss)
gen Mom_Young = (momage!=. & momage<25)
gen Dad_Black = dadrace==2
gen Dad_White = dadrace==1
gen Dad_RaceOther = (dadrace==3 | dadrace==4)
gen Dad_Hispanic = dadrace==5
gen Dad_RaceMiss = dadrace==9
gen Dad_LowEd = feduc==1
gen Dad_HighEd = feduc==2
gen Dad_EdMiss = feduc==9
gen Dad_Age = dadage
gen Dad_AgeMiss = dadagebin==9
gen Dad_AnyMissing = max(Dad_AgeMiss,Dad_RaceMiss,Dad_EdMiss)
gen Dad_Young = (dadage!=. & dadage<25)

gen FirstChild = childorder==1
gen SecondChild = childorder==2
gen ThirdPlusChild = childorder==3

gen Apgar = apgar5
gen LowBirthweight = lbw
gen MediCal = ld_med
gen Preterm = preterm_37
gen Gestation = gest
gen Male=male
gen LowIncome = lowinc
ren match_prob MatchProb

gen momedbin = .
replace momedbin = 1 if meduc==1 
replace momedbin = 2 if meduc==2
replace momedbin = 3 if meduc==9 

gen dadedbin = .
replace dadedbin = 1 if feduc==1
replace dadedbin = 2 if feduc==2
replace dadedbin = 3 if feduc==9 

gen momage_y = .
replace momage_y = 1 if momage<25 & momage!=.
replace momage_y = 2 if momage>=25 & momage!=.
replace momage_y = 3 if momage==.

gen dadage_y = .
replace dadage_y = 1 if dadage<24 & dadage!=.
replace dadage_y = 2 if dadage>=24 & dadage!=.
replace dadage_y = 3 if dadage==.

gen gestation_bin = .
replace gestation_bin = 1 if gestation_days<259 & gestation_days!=.
replace gestation_bin = 2 if gestation_days>=259 & gestation_days!=.
replace gestation_bin = 3 if gestation_days==.


gen preterm = gestation_days<259 & gestation_days!=.
gen fullterm = gestation_days>=259 & gestation_days!=.
gen gestmiss=gestation_days==.

preserve 
keep if childorder==1
*Table A9 - predicted outcomes

foreach outcome in  "numref0_2" "dayspl0_2" {

reg `outcome' i.momage_y i.momedbin i.momrace i.dadage_y i.dadedbin i.dadrace male birthweight lbw ld_med i.gestation_bin del_vagspon i.recentered_yr
predict `outcome'_hat1

}
local count = 1

local outcome  "numref0_2" 
eststo clear
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m1

	reg `outcome'_hat1 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1, vce(robust)
	sum `outcome'_hat1 if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m2
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m3
	
	reg `outcome'_hat1 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	sum `outcome'_hat1 if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m4
	
	
	esttab m1 m2 m3 m4 using "Panel`count'.tex", style(tex) replace keep(instrument) stats(ymean,labels("Outcome mean") fmt(%9.3f)) label se star(* 0.1 ** 0.05 *** 0.01) nonotes mtitles("All (Base)" "All (Predicted)" "Low Income (Base)" "Low Income (Predicted)") coeflabels(instrument "Eligible for Tax Benefits")
	
	local ++count


eststo clear
	local outcome dayspl0_2

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
	estadd scalar wmean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<`bw'
	estadd scalar vmean = r(mean)
	eststo m1
	
	reg `outcome'_hat1 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1, vce(robust)
	sum `outcome'_hat1 if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
	estadd scalar wmean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<`bw'
	estadd scalar vmean = r(mean)
	eststo m2
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
	estadd scalar wmean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<`bw'
	estadd scalar vmean = r(mean)
	eststo m3
	
	reg `outcome'_hat1 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	sum `outcome'_hat1 if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
	estadd scalar wmean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<`bw'
	estadd scalar vmean = r(mean)
	eststo m4
	
	
	esttab m1 m2 m3 m4 using "Panel8.tex", style(tex) replace keep(instrument) stats(ymean vmean wmean N,labels("Outcome mean" "Tax value mean (\\$)" "Aftertax income (\\$)" "Obs.") fmt(%9.3f 0 0 0)) label se star(* 0.1 ** 0.05 *** 0.01) nonotes mtitles("All (Base)" "All (Predicted)" "Low Income (Base)" "Low Income (Predicted)") coeflabels(instrument "Eligible for Tax Benefits")
	

panelcombine, use(Panel1.tex Panel2.tex) columncount(1) paneltitles("\# Referrals" "Days Foster Care" ) save(omnibus.tex) 

restore

*components of table A2 - Balance (manually combined)

eststo clear 

local balouts "Mom_White" "Mom_Black" "Mom_Hispanic" "Mom_RaceOther" "Mom_RaceMiss" "Mom_Young" "Mom_Age" "Mom_AgeMiss" "Mom_LowEd" "Mom_HighEd" "Mom_EdMiss" "Dad_White" "Dad_Black" "Dad_Hispanic" "Dad_RaceOther" "Dad_RaceMiss" "Dad_Young" "Dad_AgeMiss" "Dad_Age" "Dad_LowEd" "Dad_HighEd" "Dad_EdMiss" "male" "birthweight" "lbw" "ld_med" "preterm" "fullterm" "gestmiss" "del_vagspon"  "aftertax_hat" "lowinc" "MatchProb"
foreach outcome in "`balouts'" {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo `outcome'
}

esttab, se nostar stats(ymean N) keep(instrument)  

matrix C = r(coefs)
matrix S = r(stats) 
local rnames : rownames C 
local models : coleq C 
local models : list uniq models 

local i 0

foreach name of local rnames {
	local ++i 
	local j 0 
	capture matrix drop b 
	capture matrix drop se 
	foreach model of local models {
		local ++j 
		matrix tmp = C[`i', 2*`j'-1]
		if tmp[1,1]<. {
			matrix colnames tmp = `model'
			matrix b = nullmat(b), tmp
			matrix tmp[1,1] = C[`i', 2*`j']
			matrix se = nullmat(se), tmp
		}
	}
	ereturn post b
	quietly estadd matrix se 
	eststo `name'
}
local snames : rownames S
local i 0 

foreach name of local snames {
	local ++i 
	local j 0
	capture matrix drop b
	foreach model of local models {
		local ++j
		matrix tmp = S[`i', `j'] 
		matrix colnames tmp = `model'
		matrix b = nullmat(b), tmp
		
	}
	ereturn post b 
	eststo `name'
	
}
esttab instrument ymean N using "balance_allfirst.tex", style(tex) replace se mtitles("Born before Jan. 1" "Outcome mean" "Obs." ) noobs nogap compress nonumb keep("`balouts'")  star(* .1 ** .05 *** .001) coeflabels("Mom_White" "White Mom" "Mom_Black" "Black Mom" "Mom_Hispanic" "Hispanic Mom" "Mom_RaceOther" "Other Race Mom" "Mom_RaceMiss" "Mom Race Missing" "Mom_Age" "Mom's Age" "Mom_Young" "Mom $<$ 24" "Mom_AgeMiss" "Mom Age Missing" "Mom_LowEd"  "Mom Low Ed." "Mom_HighEd" "Mom High Ed." "Mom_EdMiss" "Mom Missing Ed." "Dad_White" "White Dad" "Dad_Black" "Black Dad" "Dad_Hispanic" "Hispanic Dad" "Dad_RaceOther" "Other Race Dad" "Dad_RaceMiss" "Dad Race Missing" "Dad_Age" "Dad's Age" "Dad_AgeMiss" "Dad Age Missing" "Dad_Young" "Dad $<$ 24" "Dad_LowEd" "Dad Low Ed." "Dad_HighEd" " Dad High Ed." "Dad_EdMiss" "Dad Missing Ed." "male" "Male" "birthweight" "Birthweight (g)" "lbw" "Low Birthweight" "ld_med" "Medicaid" "preterm" "Pre-term birth" "fullterm" " Full-term birth" "gestmiss" "Gestation Missing" "del_vagspon" "Spontaneous Vaginal Delivery" "aftertax_hat" "Predicted Income" "lowinc" "Low Income" "MatchProb" "Match Probability" )


eststo clear
foreach outcome in "`balouts'" {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo `outcome'
}

esttab, se nostar stats(ymean N) keep(instrument)  

matrix C = r(coefs)
matrix S = r(stats) 
local rnames : rownames C 
local models : coleq C 
local models : list uniq models 

local i 0

foreach name of local rnames {
	local ++i 
	local j 0 
	capture matrix drop b 
	capture matrix drop se 
	foreach model of local models {
		local ++j 
		matrix tmp = C[`i', 2*`j'-1]
		if tmp[1,1]<. {
			matrix colnames tmp = `model'
			matrix b = nullmat(b), tmp
			matrix tmp[1,1] = C[`i', 2*`j']
			matrix se = nullmat(se), tmp
		}
	}
	ereturn post b
	quietly estadd matrix se 
	eststo `name'
}
local snames : rownames S
local i 0 

foreach name of local snames {
	local ++i 
	local j 0
	capture matrix drop b
	foreach model of local models {
		local ++j
		matrix tmp = S[`i', `j'] 
		matrix colnames tmp = `model'
		matrix b = nullmat(b), tmp
		
	}
	ereturn post b 
	eststo `name'
	
}
esttab instrument ymean N using "balance_lowinc.tex", style(tex) nogaps replace se mtitles("Born before Jan. 1" "Outcome mean" "Obs." ) noobs nogap compress nonumb keep("`balouts'")  star(* .1 ** .05 *** .001) coeflabels("Mom_White" "White Mom" "Mom_Black" "Black Mom" "Mom_Hispanic" "Hispanic Mom" "Mom_RaceOther" "Other Race Mom" "Mom_RaceMiss" "Mom Race Missing" "Mom_Age" "Mom's Age" "Mom_Young" "Mom $<$ 24" "Mom_AgeMiss" "Mom Age Missing" "Mom_LowEd"  "Mom Low Ed." "Mom_HighEd" "Mom High Ed." "Mom_EdMiss" "Mom Missing Ed." "Dad_White" "White Dad" "Dad_Black" "Black Dad" "Dad_Hispanic" "Hispanic Dad" "Dad_RaceOther" "Other Race Dad" "Dad_RaceMiss" "Dad Race Missing" "Dad_Age" "Dad's Age" "Dad_AgeMiss" "Dad Age Missing" "Dad_Young" "Dad $<$ 24" "Dad_LowEd" "Dad Low Ed." "Dad_HighEd" " Dad High Ed." "Dad_EdMiss" "Dad Missing Ed." "male" "Male" "birthweight" "Birthweight (g)" "lbw" "Low Birthweight" "ld_med" "Medicaid"  "preterm" "Pre-term birth" "fullterm" " Full-term birth" "gestmiss" "Gestation Missing" "del_vagspon" "Spontaneous Vaginal Delivery" "aftertax_hat" "Predicted Income" "lowinc" "Low Income" "MatchProb" "Match Probability")

*Table A3 - RD estimates for all eight outcomes
*sharpened q-values calculated separately in 10_MHT.do
local count = 1
foreach outcome in  "anyref0_2"  "anyinv0_2" "anysubstal0_2" "anypl0_2" "numref0_2" "numinv0_2" "numsu0_2"  {
	eststo clear
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m1
	
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)

	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m2
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==0, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m3

	
	esttab _all using "Panel`count'.tex", style(tex) replace keep(instrument) stats(ymean,labels("Outcome mean") fmt(%9.3f)) label se star(* 0.1 ** 0.05 *** 0.01) nonotes mtitles("All" "$<$200\% FPL" "$\geq$200\% FPL") coeflabels(instrument "Eligible for Tax Benefits")
	
	local ++count
}

eststo clear
	local outcome dayspl0_2
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<`bw'
	estadd scalar vmean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
	estadd scalar wmean = r(mean)
	eststo m1
	
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
	estadd scalar wmean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<`bw'
	estadd scalar vmean = r(mean)
	eststo m2
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==0, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
	estadd scalar wmean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<`bw'
	estadd scalar vmean = r(mean)
	eststo m3
	
	esttab _all using "Panel8.tex", style(tex) replace keep(instrument) stats(ymean  vmean wmean N,labels("Outcome mean" "Tax value mean (\\$)" "Aftertax income (\\$)" "Obs.") fmt(%9.3f 0 0 0)) label se star(* 0.1 ** 0.05 *** 0.01) nonotes mtitles("All" "$<$200\% FPL" "$\geq$200\% FPL") coeflabels(instrument "Eligible for Tax Benefits")

panelcombine, use(Panel1.tex Panel2.tex Panel3.tex Panel4.tex Panel5.tex Panel6.tex Panel7.tex Panel8.tex) columncount(1) paneltitles( "Any Referrals" "Any Investigations" "Any Substantiated Referral" "Any Foster Care" "\# Referrals" "\# Investigations" "\# Substantiated Referrals" "Days Foster Care" ) save(combinedbyinc.tex) 


*Table 2 only 2 outcomes, plus adjusted p-values

local outcome "numref0_2" 
eststo clear
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1, vce(robust)
	local t = _b[instrument]/_se[instrument]
	estadd scalar pval = 2*ttail(e(df_r),abs(`t'))
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m1
	
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	local t = _b[instrument]/_se[instrument]
	estadd scalar pval = 2*ttail(e(df_r),abs(`t'))
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m2
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==0, vce(robust)
	local t = _b[instrument]/_se[instrument]
	estadd scalar pval = 2*ttail(e(df_r),abs(`t'))
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m3

	
	esttab _all using "Panel1.tex", style(tex) replace keep(instrument) stats(ymean ,labels("Outcome mean") fmt(%9.3f)) label se star(* 0.1 ** 0.05 *** 0.01) nonotes mtitles("All" "$<$200\% FPL" "$\geq$200\% FPL") coeflabels(instrument "Eligible for Tax Benefits")
	

eststo clear
	local outcome dayspl0_2
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1, vce(robust)
	local t = _b[instrument]/_se[instrument]
	estadd scalar pval = 2*ttail(e(df_r),abs(`t'))
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<`bw'
	estadd scalar vmean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
	estadd scalar wmean = r(mean)
	eststo m1
	
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	local t = _b[instrument]/_se[instrument]
	estadd scalar pval = 2*ttail(e(df_r),abs(`t'))
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
	estadd scalar wmean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<`bw'
	estadd scalar vmean = r(mean)
	eststo m2
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==0, vce(robust)
	local t = _b[instrument]/_se[instrument]
	estadd scalar pval = 2*ttail(e(df_r),abs(`t'))
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
	estadd scalar wmean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<`bw'
	estadd scalar vmean = r(mean)
	eststo m3
	
	esttab _all using "Panel2.tex", style(tex) replace keep(instrument) stats(ymean  vmean wmean N,labels("Outcome mean" "Tax value mean (\\$)" "Aftertax income (\\$)" "Obs.") fmt(%9.3f 0 0 0)) label se star(* 0.1 ** 0.05 *** 0.01) nonotes mtitles("All" "$<$200\% FPL" "$\geq$200\% FPL") coeflabels(instrument "Eligible for Tax Benefits")

panelcombine, use(Panel1.tex Panel2.tex ) columncount(1) paneltitles( "Any Referrals" "Any Investigations" "Any Substantiated Referral" "Any Foster Care" "\# Referrals" "\# Investigations" "\# Substantiated Referrals" "Days Foster Care" ) save(combinedbyincshort.tex) 

*manually add adjusted p-values 

log using "wyoungtest3000"
wyoung numref0_2 dayspl0_2,cmd(reg OUTCOMEVAR et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)) familyp(instrument) reps(3000)
log close

**Table E1 - race heterogeneity	

local outcome "numref0_2"
eststo clear
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m0
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & momrace==1 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m1

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & momrace==2 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m2
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & momrace==5 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m3

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & momrace==4 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m4
	
	
	esttab _all using "Panel1.tex", style(tex) replace keep(instrument) stats(ymean,labels("Outcome mean")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("Base" "White" "Black" "Hispanic" "Asian/PI") coeflabels(instrument "Eligible for Tax Benefits")


local outcome dayspl0_2
	
eststo clear

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
	eststo m0
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & momrace==1 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
	eststo m1

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & momrace==2 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)

sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m2


	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & momrace==5 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)

sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m3


	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & momrace==4 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)

sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m4

		
	esttab _all using "Panel2.tex", style(tex) replace keep(instrument) stats(ymean vmean wmean N,labels("Outcome mean" "Tax value mean (\\$)" "Aftertax income (\\$)" "Obs.")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("Base" "White" "Black" "Hispanic" "Asian/PI") coeflabels(instrument "Eligible for Tax Benefits")

panelcombine, use(Panel1.tex Panel2.tex) columncount(1) paneltitles("\# Referrals" "Days Foster Care") save(byrace_panelsshort.tex) 

*wyoung with subgroups - pvalues added to table manually
log using "wyoungracesubgroups"

gen newmomrace = .
replace newmomrace=momrace if momrace!=9 & momrace!=3
wyoung numref0_2 dayspl0_2,cmd(reg OUTCOMEVAR et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)) familyp(instrument) reps(3000) subgroup(newmomrace)
log close


**Table E2 - Gender heterogeneity
local outcome "numref0_2"
eststo clear
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m0
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & male==1 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m3

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & female==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m4
	
	
	esttab _all using "Panel1.tex", style(tex) replace keep(instrument) stats(ymean,labels("Outcome mean")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("Base" "Boys" "Girls") coeflabels(instrument "Eligible for Tax Benefits")


local outcome dayspl0_2
eststo clear

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
	eststo m0
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & male==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
	eststo m3

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & female==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)

sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m4	
		
	esttab _all using "Panel2.tex", style(tex) replace keep(instrument) stats(ymean vmean wmean N,labels("Outcome mean" "Tax value mean (\\$)" "Aftertax income (\\$)" "Obs.")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("Base" "Boys" "Girls") coeflabels(instrument "Eligible for Tax Benefits")

panelcombine, use(Panel1.tex Panel2.tex ) columncount(1) paneltitles("\# Referrals" "Days Foster Care") save(bygender_panelsshort.tex) 

*pvalues added to table manually
log using "wyounggendersubgroups",replace
replace male=0 if female==1
replace male=. if male==0 & female==0
preserve 
keep if indonut==1 & childorder==1 & lowinc==1
wyoung numref0_2 dayspl0_2,cmd(reg OUTCOMEVAR et et_r instrument i.recentered_yr [pw=triwt], vce(robust)) familyp(instrument) reps(3000) subgroup(male)
restore
log close


*Table A6 - poverty proxies

local outcome "numref0_2"  
eststo clear

reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m0


reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & meduc<3, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m1


reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & ld_medicaid == 1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m2
	
reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & patmiss == 1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m3
	
reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & momage<24, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m4

	esttab _all using "Panel1.tex", style(tex) replace keep(instrument) stats(ymean,labels("Outcome mean")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("Base" "Mom HS Ed." "MediCal Birth" "Father Missing Info." "Mother $<$24") coeflabels(instrument "Eligible for Tax Benefits")



eststo clear
local outcome dayspl0_2
	
reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m0

reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & meduc<3, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m1


reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & ld_medicaid == 1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m2
	
reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & patmiss == 1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m3
	
reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & momage<24, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m4
	
	esttab _all using "Panel2.tex", style(tex) replace keep(instrument) stats(ymean vmean wmean N,labels("Outcome mean" "Tax value mean (\\$)" "Aftertax income mean (\\$)" "Obs.")) label se star(* 0.1 ** 0.05 *** 0.01)  mtitles("Base" "Mom HS Ed." "MediCal Birth" "Father Missing Info." "Mother $<$24") coeflabels(instrument "Eligible for Tax Benefits")

panelcombine, use(Panel1.tex Panel2.tex) columncount(1) paneltitles("\# Referrals" "Days Foster Care") save(proxiesshort.tex) 


*Table A7 - heterogeneity by mother's country of birth
local outcome "numref0_2" 
	eststo clear
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m0
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & mom_fb==0 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m3

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & mom_fb==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m4
	
	
	esttab _all using "Panel1.tex", style(tex) replace keep(instrument) stats(ymean,labels("Outcome mean")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("Base" "Native-Born Mother" "Foreign-Born Mother") coeflabels(instrument "Eligible for Tax Benefits")

local outcome dayspl0_2
eststo clear
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
	eststo m0
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & mom_fb==0 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
	eststo m3

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & mom_fb==1 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)

sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m4	
		
	esttab _all using "Panel2.tex", style(tex) replace keep(instrument) stats(ymean vmean wmean N,labels("Outcome mean" "Tax value mean (\\$)" "Aftertax income (\\$)" "Obs.")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("Base" "Native-Born Mother" "Foreign-Born Mother") coeflabels(instrument "Eligible for Tax Benefits")

panelcombine, use(Panel1.tex Panel2.tex ) columncount(1) paneltitles("\# Referrals" "Days Foster Care") save(byfb_panelsshort.tex) 

*adjusted p-values added to table manually
log using "wyoungfbsubgroups",replace
preserve 
keep if indonut==1 & childorder==1 & lowinc==1
wyoung numref0_2 dayspl0_2,cmd(reg OUTCOMEVAR et et_r instrument i.recentered_yr [pw=triwt], vce(robust)) familyp(instrument) reps(3000) subgroup(mom_fb)
restore
log close


*Tables E3, E4 and Figures E3, E4 - allegation categories and reporter categories
gen malehirisk = male==1 & lowinc==1
gen femhirisk = female==1 & lowinc==1
gen malelorisk = male==1 & lowinc==0
gen femlorisk = female==1 & lowinc==0
gen hirisk = lowinc==1 
gen lorisk = lowinc==0
eststo clear

reg numref0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1, vce(robust)
sum numref0_2 if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
eststo repshirisk

reg numref0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & male==1, vce(robust)
sum numref0_2 if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
eststo repsmalehirisk

reg numref0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & male==0, vce(robust)
sum numref0_2 if e(sample)==1 & et<=0 & abs(et)<`bw'
estadd scalar ymean = r(mean)
eststo repsfemhirisk
	
foreach subsamp in malehirisk femhirisk hirisk {
foreach cat in physical neglect emotional other {
	reg numal`cat'0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & `subsamp'==1, vce(robust)
	sum numal`cat'0_2 if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
	estadd scalar wmean = r(mean)
	sum value if e(sample)==1 & abs(et)<`bw'
	estadd scalar vmean = r(mean)
	eststo `cat'`subsamp'
}

esttab reps`subsamp' neglect`subsamp' physical`subsamp' emotional`subsamp' other`subsamp' using "alcat`subsamp'.tex", style(tex) replace keep(instrument) stats(ymean vmean N,labels("Outcome mean" "Tax value mean (\\$)" "Obs.")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("Base" "Neglect" "Physical" "Emotional" "Other") coeflabels(instrument "Eligible for Tax Benefits")

}

panelcombine, use(alcathirisk.tex alcatmalehirisk.tex alcatfemhirisk.tex) columncount(1) paneltitles("All" "Boys" "Girls") save(alcatloedgender.tex) 

*plot 

coefplot (neglectmalehirisk,label(Boys) msymbol(S) msize(large) mlcolor(midblue) mfcolor(midblue) ciopts(lcolor(midblue))) (neglectfemhirisk,label(Girls) msymbol(D) msize(large) mlcolor(pink) mfcolor(pink) ciopts(lcolor(pink))) || (physicalmalehirisk,label(Boys)) (physicalfemhirisk,label(Girls))  || (emotionalmalehirisk,label(Boys)) (emotionalfemhirisk,label(Girls)) || (othermalehirisk,label(Boys)) (otherfemhirisk,label(Girls)) , ytitle("Effect of Eligibility for Tax Benefits",margin(medium) size(large)) xlab(1 "Neglect" 2 "Physical" 3 "Emotional" 4 "Other") vertical keep(instrument) yline(0) bycoefs scheme(s1mono)
graph export "alcatbothgenders.pdf",replace


**rep category

replace numrepother0_2 = numrepother0_2+numrepschool0_2

foreach subsamp in malehirisk femhirisk hirisk {
foreach cat in  legprof medprof cws nonmand other  {
	reg numrep`cat'0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & `subsamp'==1, vce(robust)
	sum numrep`cat'0_2 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<60
	estadd scalar wmean = r(mean)
	sum value if e(sample)==1 & abs(et)<60
	estadd scalar vmean = r(mean)
	eststo `cat'`subsamp'
}

esttab  reps`subsamp' legprof`subsamp' medprof`subsamp' cws`subsamp' nonmand`subsamp' other`subsamp' using "repcat`subsamp'.tex", style(tex) replace keep(instrument) stats(ymean vmean N,labels("Outcome mean" "Tax value mean (\\$)" "Obs.")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("Base" "Law Enforcement" "Medical Prof." "CPS caseworker" "Non-mandated" "Other") coeflabels(instrument "Eligible for Tax Benefits")

}

panelcombine, use(repcathirisk.tex repcatmalehirisk.tex repcatfemhirisk.tex) columncount(1) paneltitles("All" "Boys" "Girls") save(repcatloedgender.tex) 

*plot 

coefplot (legprofmalehirisk,label(Boys) msymbol(S) msize(large) mlcolor(midblue) mfcolor(midblue) ciopts(lcolor(midblue))) (legproffemhirisk,label(Girls) msymbol(D) msize(large) mlcolor(pink) mfcolor(pink) ciopts(lcolor(pink))) || (medprofmalehirisk,label(Boys)) (medproffemhirisk,label(Girls))  || (cwsmalehirisk,label(Boys)) (cwsfemhirisk,label(Girls)) || (nonmandmalehirisk,label(Boys)) (nonmandfemhirisk,label(Girls)) || (othermalehirisk,label(Boys)) (otherfemhirisk,label(Girls)), ytitle("Effect of Eligibility for Tax Benefits",margin(medium) size(large)) xlab(1 "Law Enfor." 2 "Medical Prof." 3 "CPS Worker" 4 "Non-Mandated" 5 "Other") vertical keep(instrument) yline(0) bycoefs scheme(s1mono)
graph export "repcatbothgenders.pdf",replace


*Figure 3a and 3b (effects by age)

foreach outcome in numref dayspl {
foreach subsamp in lorisk hirisk {
	reg `outcome'0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & `subsamp'==1 & recentered_yr<2012, vce(robust)
	eststo `outcome'`subsamp'02
	
	reg `outcome'3_5 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & `subsamp'==1 & recentered_yr<2012, vce(robust)
	eststo `outcome'`subsamp'35
	
	reg `outcome'6_8 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & `subsamp'==1 & recentered_yr<2012, vce(robust)
	eststo `outcome'`subsamp'68
}

}

foreach outcome in numref dayspl {

coefplot (`outcome'hirisk02, msize(large)) || (`outcome'hirisk35,)  || (`outcome'hirisk68,) , ytitle("Effect of Eligibility for Tax Benefits",margin(medium) size(large)) xlab(1 "Age 0 to 2" 2 "Age 3 to 5" 3 "Age 6 to 8",labsize(large)) ylab(,labsize(large)) vertical keep(instrument) yline(0) saving("agelowinc_`outcome'.gph",replace)  bycoefs scheme(s1mono) legend(off)
graph export "agelowinc_`outcome'.pdf",replace

}


*Table A5

foreach outcome in numref numinv numsu dayspl {
gen `outcome'0_8 = `outcome'0_2+`outcome'3_5+`outcome'6_8 
}
foreach outcome in ref inv  {
gen any`outcome'0_8 = (num`outcome'0_8>0 & num`outcome'0_8!=.)

}

gen anypl0_8 = (dayspl0_8>0 & dayspl0_8!=.)
gen anysu0_8 = (numsu0_8>0 & numsu0_8!=.)

local count=1
foreach outcome in  "anyref" "anyinv" "anysu" "anypl" "numref" "numinv" "numsu"    {
	eststo clear

	reg `outcome'0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered<2012, vce(robust)
	sum `outcome'0_2 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	eststo m1
	
	reg `outcome'3_5 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered<2012, vce(robust)
	sum `outcome'3_5 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	eststo m2
	
	reg `outcome'6_8 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered<2012, vce(robust)
	sum `outcome'6_8 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	eststo m3
	
	reg `outcome'0_8 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered<2012, vce(robust)
	sum `outcome'0_8 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	eststo m4
	
	
	esttab _all using "Panel`count'.tex", style(tex) replace keep(instrument) stats(ymean,labels("Outcome mean") fmt(%9.3f)) label se star(* 0.1 ** 0.05 *** 0.01) nonotes mtitles("Age 0 - 2" "Age 3 - 5" "Age 6 - 8" "Age 0 - 8") coeflabels(instrument "Eligible for Tax Benefits")
	local ++count
}

eststo clear
	local outcome dayspl
	
	reg `outcome'0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered<2012, vce(robust)
	sum `outcome'0_2 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<60
	estadd scalar vmean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<60
	estadd scalar wmean = r(mean)
	eststo m1
	
		reg `outcome'3_5 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered<2012, vce(robust)
	sum `outcome'3_5 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<60
	estadd scalar vmean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<60
	estadd scalar wmean = r(mean)
	eststo m2
	
	reg `outcome'6_8 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered<2012, vce(robust)
	sum `outcome'6_8 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<60
	estadd scalar vmean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<60
	estadd scalar wmean = r(mean)
	eststo m3
	
	reg `outcome'0_8 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered<2012, vce(robust)
	sum `outcome'0_8 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<60
	estadd scalar vmean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<60
	estadd scalar wmean = r(mean)
	eststo m4
	
	
	
	esttab _all using "Panel8.tex", style(tex) replace keep(instrument) stats(ymean vmean wmean N,labels("Outcome mean" "Tax value mean (\\$)" "Aftertax income (\\$)" "Obs.") fmt(%9.3f 0 0 0)) label se star(* 0.1 ** 0.05 *** 0.01) nonotes mtitles("Age 0 - 2" "Age 3 - 5" "Age 6 - 8" "Age 0 - 8") coeflabels(instrument "Eligible for Tax Benefits")

panelcombine, use(Panel1.tex Panel2.tex Panel3.tex Panel4.tex Panel5.tex Panel6.tex Panel7.tex Panel8.tex) columncount(1) paneltitles( "Any Referrals" "Any Investigations" "Any Substantiated Referral" "Any Foster Care" "\# Referrals" "\# Investigations" "\# Substantiated Referrals" "Days Foster Care" ) save(combinedthru8.tex) 


*Table A8

gen firstref0_2 = anyref0_2
gen firstref3_5 = (anyref0_2==0 & anyref3_5==1)
gen firstref6_8 = (anyref0_2==0 & anyref3_5==0 & anyref6_8==1) 

local count = 1
local outcome firstref
eststo clear

	reg `outcome'0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered<2012, vce(robust)
	sum `outcome'0_2 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<60
	estadd scalar vmean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<60
	estadd scalar wmean = r(mean)
	eststo m1	
	reg `outcome'3_5 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered<2012, vce(robust)
	sum `outcome'3_5 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<60
	estadd scalar vmean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<60
	estadd scalar wmean = r(mean)
	eststo m2
	
	reg `outcome'6_8 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered<2012, vce(robust)
	sum `outcome'6_8 if e(sample)==1 & et<=0 & abs(et)<60
	estadd scalar ymean = r(mean)
	sum value_hat if e(sample)==1 & abs(et)<60
	estadd scalar vmean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<60
	estadd scalar wmean = r(mean)
	eststo m3
	

	esttab _all using "firstref.tex", style(tex) replace keep(instrument) stats(ymean vmean wmean N,labels("Outcome mean" "Tax value mean (\\$)" "Aftertax income (\\$)" "Obs.") fmt(%9.3f 0 0 0)) label se star(* 0.1 ** 0.05 *** 0.01) nonotes mtitles("Age 0 - 2" "Age 3 - 5" "Age 6 - 8") coeflabels(instrument "Eligible for Tax Benefits")



*Figure A7a, A7b
eststo clear

replace hirisk = meduc==1
replace lorisk = meduc==2
foreach outcome in numref dayspl {

foreach subsamp in lorisk hirisk {
	reg `outcome'0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & `subsamp'==1 , vce(robust)
	eststo `outcome'`subsamp'02
	
	reg `outcome'0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==2 & `subsamp'==1 , vce(robust)
	eststo `outcome'`subsamp'35
	
	reg `outcome'0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder>=3 & `subsamp'==1 , vce(robust)
	eststo `outcome'`subsamp'68
	
	
}

}

coefplot (numrefhirisk02,label(Mom Low Ed.) msize(large)) (numreflorisk02,label(Mom High Ed.) msize(large)) || (numrefhirisk35,label(Mom Low Ed.)) (numreflorisk35,label(Mom High Ed.))  || (numrefhirisk68,label(Mom Low Ed.)) (numreflorisk68,label(Mom High Ed.)) , ytitle("Effect of Eligibility for Tax Benefits",margin(medium) size(large)) xlab(1 "1st Child" 2 "2nd Child" 3 "3rd+ Child",labsize(large)) ylab(,labsize(large)) vertical keep(instrument) yline(0) saving("childorderlohi_numref.gph",replace)  bycoefs scheme(s1mono) legend(off)
graph export "childorderlohi_numref.pdf",replace


coefplot (daysplhirisk02,label(Mom Low Ed.) msize(large)) (dayspllorisk02,label(Mom High Ed.) msize(large)) || (daysplhirisk35,label(Mom Low Ed.)) (dayspllorisk35,label(Mom High Ed.))  || (daysplhirisk68,label(Mom Low Ed.)) (dayspllorisk68,label(Mom High Ed.)) , ytitle("Effect of Eligibility for Tax Benefits",margin(medium) size(large)) xlab(1 "1st Child" 2 "2nd Child" 3 "3rd+ Child",labsize(large)) ylab(,labsize(large)) vertical keep(instrument) yline(0) saving("childorderlohi_dayspl.gph",replace) bycoefs scheme(s1mono) legend(size(large))
graph export "childorderlohi_dayspl.pdf",replace


*Table C2 - effects before vs. after paid family leave 
gen prepfl = recentered_yr<=2004
local count = 1
local outcome "numref0_2" 
	eststo clear
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m0
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & prepfl==1 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m1

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & prepfl==0 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	eststo m2
	
	esttab _all using "Panel`count'.tex", style(tex) replace keep(instrument) stats(ymean,labels("Outcome mean")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("Base" "Pre PFL" "Post PFL") coeflabels(instrument "Eligible for Tax Benefits")
	local ++count


local outcome dayspl0_2
	
eststo clear

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
	eststo m0
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & prepfl==1 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)
	sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
	eststo m1

	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & prepfl==0 & lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)

sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m2

		
	esttab _all using "Panel2.tex", style(tex) replace keep(instrument) stats(ymean vmean wmean N,labels("Outcome mean" "Tax value mean (\\$)" "Aftertax income (\\$)" "Obs.")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("Base" "Pre PFL" "Post PFL") coeflabels(instrument "Eligible for Tax Benefits")

panelcombine, use(Panel1.tex Panel2.tex) columncount(1) paneltitles("\# Referrals" "Days Foster Care") save(het_pfl.tex) 

***Table A4 - distribution of foster care days 

gen dayspl_0 = dayspl0_2>0
gen dayspl_7 = dayspl0_2>7
gen dayspl_30 = dayspl0_2>30
gen dayspl_180 = dayspl0_2>180
gen dayspl_360=dayspl0_2>360

eststo clear
local i = 1

foreach outcome in  "dayspl_0"  "dayspl_7" "dayspl_30" "dayspl_180" "dayspl_360"  {
	
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 &  lowinc==1, vce(robust)
	sum `outcome' if e(sample)==1 & et<=0 & abs(et)<`bw'
	estadd scalar ymean = r(mean)

sum aftertax if e(sample)==1 & abs(et)<`bw'
estadd scalar wmean = r(mean)
sum value if e(sample)==1 & abs(et)<`bw'
estadd scalar vmean = r(mean)
eststo m`i'

	local ++i
}

	esttab _all using "dayspl_distribution.tex", style(tex) replace keep(instrument) stats(ymean vmean wmean N,labels("Outcome mean" "Tax value mean (\\$)" "Aftertax income (\\$)" "Obs.")) label se star(* 0.1 ** 0.05 *** 0.01) mtitles("$>$0" "$>$7" "$>$30" "$>$180" "$>$360") coeflabels(instrument "Eligible for Tax Benefits")
	
***binscatters***	

local tick 10
replace et = et*-1
local outcome Mom_Black 
	local yl = 0
	local yu = 0.1
	local yt = 0.01
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	

	local yl = 0.03
	local yu = 0.13
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	

local outcome Mom_Age 
	local yl = 25
	local yu = 27
	local yt = .2
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	

	local yl = 21
	local yu = 23

binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	

local outcome Mom_LowEd 
	local yl = 0.39
	local yu = 0.49
	local yt = 0.01
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	
	local yl = 0.7
	local yu = 0.8
	
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	

local outcome Mom_AnyMissing
	local yl = 0.02
	local yu = 0.12
	local yt = 0.01
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	

binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	

local outcome Dad_Black 
	local yl = 0
	local yu = 0.1
	local yt = 0.01
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	

binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	

local outcome Dad_Age 
	local yl = 28
	local yu = 30
	local yt = .2
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	
	local yl = 23.4
	local yu = 25.4
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	

local outcome Dad_LowEd 
	local yl = 0.39
	local yu = 0.49
	local yt = 0.01
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	

	local yl = 0.59
	local yu = 0.69
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	

local outcome Dad_AnyMissing
	local yl = 0.12
	local yu = 0.22
	local yt = 0.01
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	
	local yl = 0.25
	local yu = 0.35
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	


local outcome male
	local yl = 0.45
	local yu = 0.55
	local yt = 0.01
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	

binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	

local outcome birthweight
	local yl = 3100
	local yu = 3400
	local yt = 50
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	

binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	


local outcome ld_med
	local yl = 0.35
	local yu = 0.45
	local yt = 0.01
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	
	local yl = 0.64
	local yu = 0.74
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	

local outcome preterm
	local yl = 0.05
	local yu = 0.15
	local yt = 0.01
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	

binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	

local outcome del_vagspon
	local yl = 0.6
	local yu = 0.75
	local yt = 0.02
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	

binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	

local outcome gestation
	local yl = 270
	local yu = 280
	local yt = 1
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	

binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	


local outcome aftertax_hat
	local yl = 56000
	local yu = 66000
	local yt = 1000
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	
	local yl = 20000
	local yu = 28000
	local yt = 1000
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	


local outcome MatchProb
	local yl = 0.98
	local yu = 1
	local yt = 0.005
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'.pdf", replace	

binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "val`outcome'lowinc.pdf", replace	


*first stage 
gen firststage = aftertax_hat
replace firststage = aftertax_hat+value_hat if et<0

local outcome firststage 
local yl = 58000
	local yu = 68000
	local yt = 1000
	*local ytitle = "Num. referrals"
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("Predicted After-tax Income")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "binscatter`outcome'.pdf", replace

local yl = 20000
	local yu = 30000
	local yt = 1000
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("Predicted After-tax Income")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "binscatterlowinc`outcome'.pdf", replace


local outcome numref0_2 
	local yl = .15
	local yu = 0.35
	local yt = 0.05
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) color(black black) legend(off)
graph export "binscatterlowinc`outcome'.pdf", replace


binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols( circle  square) color(black gray) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) by(male) lcolors(black gray) legend(order(1 "Female" 2 "Male"))
graph export "binscatterlowinc`outcome'gender.pdf", replace

	local yl = .05
	local yu = 0.45
	local yt = 0.05
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle square) color(black  gray) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) by(mom_fb) lcolors( black  gray) legend(order(1 "Native-born mother" 2 "Foreign-born mother"))
graph export "binscatterlowinc`outcome'fb.pdf", replace


gen racebw =.
replace racebw = 1 if  momrace==1
replace racebw = 2 if  momrace==2 


	local yl = .2
	local yu = 0.8
	local yt = 0.1
	
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols( circle  square) color(black gray) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) by(racebw) lcolors( black  gray) legend(order(1 "White mother" 2 "Black mother"))
graph export "binscatterlowinc`outcome'race.pdf", replace

local outcome numinv0_2 
	local yl = .1
	local yu = 0.3
	local yt = 0.05
binscatter `outcome' et if abs(et)<`bw' & childorder==1 &  lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge)) ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) color(black black)  ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large))  lcolors(black) legend(off)
graph export "binscatterlowinc`outcome'.pdf", replace

local outcome numsu0_2 
	local yl = 0.02
	local yu = 0.08
	local yt = 0.01
binscatter `outcome' et if abs(et)<`bw' & childorder==1 &  lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge)) ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) color(black black)  ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large))  lcolors(black) legend(off)
graph export "binscatterlowinc`outcome'.pdf", replace


local outcome dayspl0_2 
	local yl = 2
	local yu = 16
	local yt = 2
binscatter `outcome' et if abs(et)<`bw' & childorder==1 &  lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge)) ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) color(black black)  ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large))  lcolors(black) legend(off)
graph export "binscatterlowinc`outcome'.pdf", replace


binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols( circle  square) color(black gray) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) by(male) lcolors(black gray) legend(order(1 "Female" 2 "Male"))
graph export "binscatterlowinc`outcome'gender.pdf", replace

	local yl = 0
	local yu = 16
	local yt = 2
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle square) color(black  gray) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) by(mom_fb) lcolors( black  gray) legend(order(1 "Native-born mother" 2 "Foreign-born mother"))
graph export "binscatterlowinc`outcome'fb.pdf", replace


	local yl = 0
	local yu = 24
	local yt = 4

binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols( circle  square) color(black gray) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) by(racebw) lcolors( black  gray) legend(order(1 "White mother" 2 "Black mother"))
graph export "binscatterlowinc`outcome'race.pdf", replace

local outcome anyref0_2 
	local yl = .08
	local yu = 0.2
	local yt = 0.04
binscatter `outcome' et if abs(et)<`bw' & childorder==1 &  lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge)) ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) color(black black)  ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large))  lcolors(black) legend(off)
graph export "binscatterlowinc`outcome'.pdf", replace

local outcome anyinv0_2 
	local yl = .08
	local yu = 0.2
	local yt = 0.04
binscatter `outcome' et if abs(et)<`bw' & childorder==1 &  lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge)) ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) color(black black)  ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large))  lcolors(black) legend(off)
graph export "binscatterlowinc`outcome'.pdf", replace


local outcome anysu0_2 
	local yl = 0.02
	local yu = .08
	local yt = .02
binscatter `outcome' et if abs(et)<`bw' & childorder==1 &  lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge)) ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) color(black black)  ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large))  lcolors(black) legend(off)
graph export "binscatterlowinc`outcome'.pdf", replace


local outcome anypl0_2 
	local yl = .0
	local yu = .04
	local yt = .01
binscatter `outcome' et if abs(et)<`bw' & childorder==1 &  lowinc==1 & indonut==1, discrete line(lfit) xline(0) xlabel(`lower'(`tick')`upper',labsize(medlarge)) ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) color(black black)  ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large))  lcolors(black) legend(off)
graph export "binscatterlowinc`outcome'.pdf", replace

***0 to 8 results
foreach var in "numref" "dayspl" {
gen `var'0_8 = `var'0_2+`var'3_5+`var'6_8
}

local outcome numref0_8 
	local yl = .55
	local yu = 0.85
	local yt = 0.05
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & recentered<2012 & indonut==1,  discrete xline(0) line(lfit) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) color(black) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) lcolors( black) legend(off)
graph export "binscatterlowinc`outcome'.pdf", replace

local outcome dayspl0_8
	local yl = 10
	local yu = 30
	local yt = 5
binscatter `outcome' et if abs(et)<`bw' & childorder==1 & lowinc==1 & recentered<2012 & indonut==1,  discrete xline(0) line(lfit) xlabel(`lower'(`tick')`upper',labsize(medlarge))  ylabel(`yl'(`yt')`yu',labsize(medlarge)) yscale(r(`yl' `yu')) rd(0) msymbols(circle) color(black) ytitle("`ytitle'")  xtitle("Days to Jan 1",size(large)) lcolors( black) legend(off)
graph export "binscatterlowinc`outcome'.pdf", replace


*Figure 1 - number of births by day
preserve
gen num_first = (childorder==1)
collapse (sum) num_first,by(et birth_y indonut)
drop if birth_y==.
collapse (mean) numbirths num_first num_1stlowinc,by(et indonut)


local var num_first
local yl 0
local yu 750
local yt 250

twoway (bar `var' et if indonut==0 & abs(et)<=60, color(gs12)) (bar `var' et if indonut==1 & abs(et)<=60, xline(0) ylabel(`yl'(`yt')`yu') yscale(r(`yl' `yu')) xlabel(`lower'(`tick')`upper') graphregion(color(white)) bgcolor(white) color( gs4) legend(off) xtitle("Days to Jan 1") ytitle("Average Number of Births"))
graph export "`var'_donut.pdf", as(pdf) replace


***Robustness
*Figure A10a,b (by excluded year)
replace et = et*-1

local i = 1
foreach outcome in numref dayspl  {

forval yr = 2000/2017 {
	reg `outcome'0_2 et et_r instrument [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & recentered_yr!=`yr'
	eststo m`yr'
}


if `i'==1 {
	local yline -0.0206

}

if `i'==2 {
	local yline -1.666

}

coefplot m2000 || m2001 || m2002 || m2003 || m2004 || m2005 || m2006 || m2007 || m2008 || m2009 || m2010 || m2011 || m2012 || m2013 || m2014 || m2015 || m2016 || m2017, vertical keep(instrument) yline(0,lcolor(black)) yline(`yline', lcolor(red) lpattern(dash)) xtitle("Recentered Year") ytitle("Effect of Eligibility for Tax Benefits") msymbol(circle_hollow) msize(large) mcolor(black) color(black) bycoefs xlab(1 "2000" 2 "2001" 3 "2002" 4 "2003" 5 "2004" 6 "2005" 7 "2006" 8 "2007" 9 "2008" 10 "2009" 11 "2010" 12 "2011" 13 "2012" 14 "2013" 15 "2014" 16 "2015" 17 "2016" 18 "2017", labsize(medsmall) angle(45)) scheme(s1mono)
graph export "robbyyearex`outcome'.pdf", as(pdf) replace
local ++i
}

*Figure 9a,b - bandwidth and donut robustness 
*optimal bandwidth from 11_Bandwidth.do
eststo clear
local i 1
foreach outcome in numref dayspl {

forval bw = 30(10)90 {
local lower -`bw'
local upper `bw'
replace transformed = et/(`bw')
replace inbandwidth = abs(et)<`bw'
replace triwt = (1-abs(transformed))*inbandwidth

foreach donut of numlist 6 8 10 12 14 {
	
	local donutl `donut'
	local donutr `donut' 
	
	
	replace indonut = 0
	replace indonut = 1 if et<-(`donutl'-1) | et>(`donutr')

	if `bw'> `donut' {
		reg `outcome'0_2 et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 , vce(robust)
		eststo bw`bw'd`donut'
	}
	}
}

if `i'==1 {
	local yline -0.0206
	local opt = -.02346

}

if `i'==2 {
	local yline -1.666
	local opt = -2.82403

}

coefplot (bw30d6,label(6) msymbol(S) mlcolor(midgreen) mfcolor(midgreen) ciopts(lcolor(midgreen))) (bw30d8,label(8) msymbol(X) mlcolor(black) mfcolor(black) ciopts(lcolor(black))) (bw30d10,label(10) msymbol(T) mlcolor(purple) mfcolor(purple) ciopts(lcolor(purple))) (bw30d12,label(12) msymbol(O) mlcolor(red) mfcolor(red) ciopts(lcolor(red))) (bw30d14,label(14) msymbol(D) mlcolor(midblue) mfcolor(midblue) ciopts(lcolor(midblue)))    || (bw40d6,label(6) msymbol(S) mlcolor(midgreen) mfcolor(midgreen) ciopts(lcolor(midgreen))) (bw40d8,label(8) msymbol(X) mlcolor(black) mfcolor(black) ciopts(lcolor(black))) (bw40d10,label(10) msymbol(T) mlcolor(purple) mfcolor(purple) ciopts(lcolor(purple))) (bw40d12,label(12) msymbol(O) mlcolor(red) mfcolor(red) ciopts(lcolor(red))) (bw40d14,label(14) msymbol(D) mlcolor(midblue) mfcolor(midblue) ciopts(lcolor(midblue))) || (bw50d6,label(6) msymbol(S) mlcolor(midgreen) mfcolor(midgreen) ciopts(lcolor(midgreen))) (bw50d8,label(8) msymbol(X) mlcolor(black) mfcolor(black) ciopts(lcolor(black))) (bw50d10,label(10) msymbol(T) mlcolor(purple) mfcolor(purple) ciopts(lcolor(purple))) (bw50d12,label(12) msymbol(O) mlcolor(red) mfcolor(red) ciopts(lcolor(red))) (bw50d14,label(14) msymbol(D) mlcolor(midblue) mfcolor(midblue) ciopts(lcolor(midblue))) || (bw60d6,label(6) msymbol(S) mlcolor(midgreen) mfcolor(midgreen) ciopts(lcolor(midgreen))) (bw60d8,label(8) msymbol(X) mlcolor(black) mfcolor(black) ciopts(lcolor(black))) (bw60d10,label(10) msymbol(T) mlcolor(purple) mfcolor(purple) ciopts(lcolor(purple))) (bw60d12,label(12) msymbol(O) mlcolor(red) mfcolor(red) ciopts(lcolor(red))) (bw60d14,label(14) msymbol(D) mlcolor(midblue) mfcolor(midblue) ciopts(lcolor(midblue))) || (bw70d6,label(6) msymbol(S) mlcolor(midgreen) mfcolor(midgreen) ciopts(lcolor(midgreen))) (bw70d8,label(8) msymbol(X) mlcolor(black) mfcolor(black) ciopts(lcolor(black))) (bw70d10,label(10) msymbol(T) mlcolor(purple) mfcolor(purple) ciopts(lcolor(purple))) (bw70d12,label(12) msymbol(O) mlcolor(red) mfcolor(red) ciopts(lcolor(red))) (bw70d14,label(14) msymbol(D) mlcolor(midblue) mfcolor(midblue) ciopts(lcolor(midblue))) || (bw80d6,label(6) msymbol(S) mlcolor(midgreen) mfcolor(midgreen) ciopts(lcolor(midgreen))) (bw80d8,label(8) msymbol(X) mlcolor(black) mfcolor(black) ciopts(lcolor(black))) (bw80d10,label(10) msymbol(T) mlcolor(purple) mfcolor(purple) ciopts(lcolor(purple))) (bw80d12,label(12) msymbol(O) mlcolor(red) mfcolor(red) ciopts(lcolor(red))) (bw80d14,label(14) msymbol(D) mlcolor(midblue) mfcolor(midblue) ciopts(lcolor(midblue))) || (bw90d6,label(6) msymbol(S) mlcolor(midgreen) mfcolor(midgreen) ciopts(lcolor(midgreen))) (bw90d8,label(8) msymbol(X) mlcolor(black) mfcolor(black) ciopts(lcolor(black))) (bw90d10,label(10) msymbol(T) mlcolor(purple) mfcolor(purple) ciopts(lcolor(purple))) (bw90d12,label(12) msymbol(O) mlcolor(red) mfcolor(red) ciopts(lcolor(red))) (bw90d14,label(14) msymbol(D) mlcolor(midblue) mfcolor(midblue) ciopts(lcolor(midblue))),  xlab(1 "30" 2 "40" 3 "50" 4 "60" 5 "70" 6 "80" 7 "90") vertical keep(instrument) yline(0) yline(`yline',lpattern(dash)) yline(`opt',lpattern(longdash))  bycoefs xtitle("Bandwidth Size (days)",margin(medium) size(large)) scheme(s1mono) legend(rows(1))
graph export "bwd`outcome'.pdf", as(pdf) replace


*Figure A11a,b (other holidays)
eststo clear

foreach outcome in  "numref0_2" "dayspl0_2"   {
gen b`outcome'=.
gen s`outcome'=.
}

gen holidaylist = ""
replace holidaylist = "January 1" if _n==1
replace holidaylist = "MLK Day" if _n==2
replace holidaylist = "President's Day" if _n==3
replace holidaylist = "Memorial Day" if _n==4
replace holidaylist = "July 4" if _n==5
replace holidaylist = "Labor Day" if _n==6
replace holidaylist = "Columbus Day" if _n==7
replace holidaylist = "Veterans Day" if _n==8
replace holidaylist = "Thanksgiving" if _n==9
replace holidaylist = "Christmas Day" if _n==10


foreach outcome in "numref0_2" "dayspl0_2"   {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 
	replace b`outcome' = _b[instrument] if holidaylist=="January 1"
	replace s`outcome' = _se[instrument] if holidaylist=="January 1"
}


gen dow = dow(bdate_mdy)
gen xmas = mdy(12,25,birth_y)
gen ny = mdy(1,1,birth_y)
replace ny = mdy(1,1,birth_y+1) if birth_m==12
gen j4 = mdy(7,4,birth_y)
gen vet = mdy(11,11,birth_y)
gen ontg = .
replace ontg=bdate_mdy if dow==4 & birth_m==11 & birth_d>=22 & birth_d<=28
bys birth_y:egen tg=mode(ontg)
tab tg
tab birth_y if tg==.
gen onme = .
replace onme=bdate_mdy if dow==1 & birth_m==5 & birth_d>=25 & birth_d<=31
bys birth_y:egen mem=mode(onme)
gen onla = .
replace onla=bdate_mdy if dow==1 & birth_m==9 & birth_d>=1 & birth_d<=7
bys birth_y:egen lab=mode(onla)

gen onmlk = .
replace onmlk=bdate_mdy if dow==1 & birth_m==1 & birth_d>=15 & birth_d<=21
bys birth_y:egen mlk=mode(onmlk)

gen onpr = .
replace onpr=bdate_mdy if dow==1 & birth_m==2 & birth_d>=15 & birth_d<=21
bys birth_y:egen pre=mode(onpr)

gen oncol = .
replace oncol=bdate_mdy if dow==1 & birth_m==10 & birth_d>=8 & birth_d<=14
bys birth_y:egen col=mode(oncol)


drop ontg onme onla onmlk onpr oncol

*tg
forval x = 1999/2017 {
	gen tg`x' = .
	replace tg`x' = tg if birth_y==`x'
}
forval x = 1999/2017 {
	egen tgyr`x' = max(tg`x')
}

forval x = 1999/2017 {
	drop tg`x'
}

gen tgdate = tg

forval x = 2000/2017 {
	local y = `x'-1
	replace tgdate = tgyr`y' if birth_m<7 & birth_y==`x'
}

replace recentered_yr = year(tgdate)
replace et = (bdate_mdy - tgdate)*-1
replace instrument = (et>=0)
replace et_r = et*instrument
replace transformed = et/(`bw')
replace inbandwidth = abs(et)<`bw'
replace triwt = (1-abs(transformed))*inbandwidth
replace indonut = 0
replace indonut = 1 if et<-(`donutl'-1) | et>(`donutr')

foreach outcome in   "numref0_2" "dayspl0_2"   {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 
	replace b`outcome' = _b[instrument] if holidaylist=="Thanksgiving"
	replace s`outcome' = _se[instrument] if holidaylist=="Thanksgiving"
}


*j4 
gen jul4date = mdy(7,4,birth_y)
replace recentered_yr = year(jul4date)

replace et = (bdate_mdy - jul4date)*-1
replace instrument = (et>=0)
replace et_r = et*instrument
replace transformed = et/(`bw')
replace inbandwidth = abs(et)<`bw'
replace triwt = (1-abs(transformed))*inbandwidth
replace indonut = 0
replace indonut = 1 if et<-(`donutl'-1) | et>(`donutr')

foreach outcome in   "numref0_2" "dayspl0_2"   {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & birth_y>1999 & birth_y<2017
	replace b`outcome' = _b[instrument] if holidaylist=="July 4"
	replace s`outcome' = _se[instrument] if holidaylist=="July 4"
}


*memorial 
replace recentered_yr = year(mem)

replace et = (bdate_mdy - mem)*-1
replace instrument = (et>=0)
replace et_r = et*instrument
replace transformed = et/(`bw')
replace inbandwidth = abs(et)<`bw'
replace triwt = (1-abs(transformed))*inbandwidth
replace indonut = 0
replace indonut = 1 if et<-(`donutl'-1) | et>(`donutr')

foreach outcome in   "numref0_2" "dayspl0_2"   {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 & birth_y>1999 & birth_y<2017
	replace b`outcome' = _b[instrument] if holidaylist=="Memorial Day"
	replace s`outcome' = _se[instrument] if holidaylist=="Memorial Day"
}


*labor 
replace recentered_yr = year(lab)

replace et = (bdate_mdy - lab)*-1
replace instrument = (et>=0)
replace et_r = et*instrument
replace transformed = et/(`bw')
replace inbandwidth = abs(et)<`bw'
replace triwt = (1-abs(transformed))*inbandwidth
replace indonut = 0
replace indonut = 1 if et<-(`donutl'-1) | et>(`donutr')

foreach outcome in   "numref0_2" "dayspl0_2"   {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 
	replace b`outcome' = _b[instrument] if holidaylist=="Labor Day"
	replace s`outcome' = _se[instrument] if holidaylist=="Labor Day"
}



*columbus 
replace recentered_yr = year(col)

replace et = (bdate_mdy - col)*-1
replace instrument = (et>=0)
replace et_r = et*instrument
replace transformed = et/(`bw')
replace inbandwidth = abs(et)<`bw'
replace triwt = (1-abs(transformed))*inbandwidth
replace indonut = 0
replace indonut = 1 if et<-(`donutl'-1) | et>(`donutr')

foreach outcome in   "numref0_2" "dayspl0_2"   {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 
	replace b`outcome' = _b[instrument] if  holidaylist=="Columbus Day"
	replace s`outcome' = _se[instrument] if holidaylist=="Columbus Day"
}


*veterans
gen vetdate = vet
replace vetdate = mdy(11,11,birth_y+1) if birth_m<6

replace recentered_yr = year(vetdate)

replace et = (bdate_mdy - vetdate)*-1
replace instrument = (et>=0)
replace et_r = et*instrument
replace transformed = et/(`bw')
replace inbandwidth = abs(et)<`bw'
replace triwt = (1-abs(transformed))*inbandwidth
replace indonut = 0
replace indonut = 1 if et<-(`donutl'-1) | et>(`donutr')

foreach outcome in   "numref0_2" "dayspl0_2"   {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 
	replace b`outcome' = _b[instrument] if  holidaylist=="Veterans Day"
	replace s`outcome' = _se[instrument] if holidaylist=="Veterans Day"
}



*christmas
gen xmdate = mdy(12,25,birth_y)
replace xmdate = mdy(12,25,birth_y-1) if birth_m<6

replace recentered_yr = year(xmdate)

replace et = (bdate_mdy - xmdate)*-1
replace instrument = (et>=0)
replace et_r = et*instrument
replace transformed = et/(`bw')
replace inbandwidth = abs(et)<`bw'
replace triwt = (1-abs(transformed))*inbandwidth
replace indonut = 0
replace indonut = 1 if et<-(`donutl'-1) | et>(`donutr')

foreach outcome in   "numref0_2" "dayspl0_2"   {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 
	replace b`outcome' = _b[instrument] if  holidaylist=="Christmas Day"
	replace s`outcome' = _se[instrument] if holidaylist=="Christmas Day"
}

*mlk 
forval x = 1999/2017 {
	gen mlk`x' = .
	replace mlk`x' = mlk if birth_y==`x'
}
forval x = 1999/2017 {
	egen mlkyr`x' = max(mlk`x')
}

forval x = 1999/2017 {
	drop mlk`x'
}

gen mlkdate = mlk

forval x = 1999/2016 {
	local y = `x'+1
	replace mlkdate = mlkyr`y' if birth_m>6 & birth_y==`x'
}

replace recentered_yr = year(mlkdate)
replace et = (bdate_mdy - mlkdate)*-1
replace instrument = (et>=0)
replace et_r = et*instrument
replace transformed = et/(`bw')
replace inbandwidth = abs(et)<`bw'
replace triwt = (1-abs(transformed))*inbandwidth
replace indonut = 0
replace indonut = 1 if et<-(`donutl'-1) | et>(`donutr')
tab recentered
foreach outcome in   "numref0_2" "dayspl0_2"   {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 
	replace b`outcome' = _b[instrument] if  holidaylist=="MLK Day"
	replace s`outcome' = _se[instrument] if holidaylist=="MLK Day"
}


*pres 
forval x = 1999/2017 {
	gen pres`x' = .
	replace pres`x' = pre if birth_y==`x'
}
forval x = 1999/2017 {
	egen presyr`x' = max(pres`x')
}

forval x = 1999/2017 {
	drop pres`x'
}

gen presdate = pre

forval x = 1999/2016 {
	local y = `x'+1
	replace presdate = presyr`y' if birth_m>6 & birth_y==`x'
}

replace recentered_yr = year(presdate)

replace et = (bdate_mdy - presdate)*-1
replace instrument = (et>=0)
replace et_r = et*instrument
replace transformed = et/(`bw')
replace inbandwidth = abs(et)<`bw'
replace triwt = (1-abs(transformed))*inbandwidth
replace indonut = 0
replace indonut = 1 if et<-(`donutl'-1) | et>(`donutr')
foreach outcome in   "numref0_2" "dayspl0_2"   {
	reg `outcome' et et_r instrument i.recentered_yr [pw=triwt] if indonut==1 & childorder==1 & lowinc==1 
	replace b`outcome' = _b[instrument] if  holidaylist=="President's Day"
	replace s`outcome' = _se[instrument] if holidaylist=="President's Day"
}


preserve
keep banyref0_2 banyinv0_2 banysu0_2 banypl0_2 bnumref0_2 bnuminv0_2 bnumsu0_2 bdayspl0_2 sanyref0_2 sanyinv0_2 sanysu0_2 sanypl0_2 snumref0_2 snuminv0_2 snumsu0_2 sdayspl0_2 holidaylist 
duplicates drop 
drop if holiday==""
gen num=_n
foreach outcome in   "numref0_2" "dayspl0_2"   {
gen cl_`outcome' = b`outcome'-1.96*s`outcome' 
gen cu_`outcome' = b`outcome'+1.96*s`outcome'
}
save "holidayestimate.dta",replace 
use "holidayestimate.dta",clear 

sort num

foreach outcome in   "numref0_2" "dayspl0_2"   {
	local jan1beta = b`outcome'[1]
	local ytext = -2*`jan1beta'
twoway (scatter b`outcome' num) (rcap cl_`outcome' cu_`outcome' num, yline(`jan1beta',lcolor(red)) yline(0,lcolor(black)) xtitle("") xlab(1 "January 1" 2 "MLK Day" 3 "President's Day" 4 "Memorial Day" 5 "July 4" 6 "Labor Day" 7 "Columbus Day" 8 "Veteran's Day" 9 "Thanksgiving" 10 "Christmas",angle(45)) scale(1.5) scheme(s1mono) legend(off) )
graph export "holiday placebo `outcome'.pdf", as(pdf) name("Graph") replace

	}

restore
